# Rossler模型构造
import matplotlib.pyplot as plt
import numpy as np
from pylab import mpl
"""
# 罗斯勒方程
dx/dt = - y - z
dy/dt = x + a * y
dz/dt = b + z(x - c)
"""

# 模型参数
a = 0.5
b = 2.0
c = 4.0
t = [0]
x = [1]
y = [0]
z = [0]
h = 0.01
n = 20000

mpl.rcParams['font.sans-serif'] = ['FangSong'] # 指定默认字体
mpl.rcParams['axes.unicode_minus'] = False # 解决保存图像是负号'-'显示为方块的问题

def Rossler(x0, y0, z0, w, p, q, r, T):
    h = 0.01
    x = []
    y = []
    z = []
    for t in range(T):
        xt = x0 + h * (- w * y0 - z0)
        yt = y0 + h * (w * x0 + p * y0)
        zt = z0 + h * (q + z0 * (x0 - r))

        #x0、y0、z0统一更新
        x0, y0, z0 = xt, yt, zt
        x.append(x0)
        y.append(y0)
        z.append(z0)

    return x, y, z


# 辅助龙格库塔方程组进行迭代构造全部数据
# dx/dt = - y - z
def FX(X, Y, Z):
    return - Y - Z

# dy/dt = x + a * y
def FY(X, Y, Z):
    return X + a * Y

# dz/dt = b + z(x - c)
def FZ(X, Y, Z):
    return b + Z * (X - c)

# 利用龙格库塔方程组对Rossler微方程组求解
def main():
    for i in range(n):
        K1 = FX(x[-1], y[-1], z[-1])
        L1 = FY(x[-1], y[-1], z[-1])
        M1 = FZ(x[-1], y[-1], z[-1])
        K2 = FX(x[-1] + h * K1 / 2, y[-1] + h * L1 / 2, z[-1] + h * M1 / 2)
        L2 = FY(x[-1] + h * K1 / 2, y[-1] + h * L1 / 2, z[-1] + h * M1 / 2)
        M2 = FZ(x[-1] + h * K1 / 2, y[-1] + h * L1 / 2, z[-1] + h * M1 / 2)
        K3 = FX(x[-1] + h * K2 / 2, y[-1] + h * L2 / 2, z[-1] + h * M2 / 2)
        L3 = FY(x[-1] + h * K2 / 2, y[-1] + h * L2 / 2, z[-1] + h * M2 / 2)
        M3 = FZ(x[-1] + h * K2 / 2, y[-1] + h * L2 / 2, z[-1] + h * M2 / 2)
        K4 = FX(x[-1] + h * K3, y[-1] + h * L3, z[-1] + h * M3)
        L4 = FY(x[-1] + h * K3, y[-1] + h * L3, z[-1] + h * M3)
        M4 = FZ(x[-1] + h * K3, y[-1] + h * L3, z[-1] + h * M3)
        x.append(x[-1] + h / 6 * (K1 + 2 * K2 + 2 * K3 + K4))
        y.append(y[-1] + h / 6 * (L1 + 2 * L2 + 2 * L3 + L4))
        z.append(z[-1] + h / 6 * (M1 + 2 * M2 + 2 * M3 + M4))


    # 构造Rossler模型三维演化轨迹
    fig = plt.figure()
    ax = fig.add_subplot(projection='3d')
    ax.plot(x, y, z, label='parametric curve', c='orange')
    plt.title('Rossler模型三维演化轨迹')
    plt.xlabel('x轴')
    plt.ylabel('y轴')
    #plt.savefig('Rossler三维演化轨迹.png')
    plt.show()

    #设置参量
    w = 1
    p = 0.165
    q = 0.2
    # 给出迭代前的初值
    x0 = 1
    y0 = 0
    z0 = 0
    x1,y1,z1 = Rossler(x0, y0, z0, w, p, q, 10, n)
    #初值微小的变化
    x2 = 1
    y2 = 0
    z2 = 0.00001
    xx,yy,zz = Rossler(x2, y2, z2, w, p, q, 10, n)
    t = np.arange(0, n)


    # x(t)轨道演化图
    #plt.scatter(t, x1, s=1, c="green")
    #plt.scatter(t, xx, s=1,)
    #plt.title('x(t)轨道演化图')
    #plt.xlabel('x轴')
    #plt.ylabel('y轴')
    #plt.savefig('Rossler的x(t)轨道演化图.png')

    # y(t)轨道演化图
    #plt.scatter(t, y1, s=1, c='pink')
    #plt.scatter(t, yy, s=1)
    #plt.title('y(t)轨道演化图')
    #plt.xlabel('x轴')
    #plt.ylabel('y轴')
    #plt.savefig('Rossler的y(t)轨道演化图.png')

    # z(t)轨道演化图
    #plt.plot(t, z1,  c='purple',linestyle='--' )
    #plt.plot(t, zz, )
    plt.scatter(t, z1, s=1, c='purple')
    plt.scatter(t, zz, s=1)
    plt.title('z(t)轨道演化图')
    plt.xlabel('x轴')
    plt.ylabel('y轴')
    plt.savefig('Rossler的z(t)轨道演化图.png')

    plt.show()


if __name__ == '__main__':
    main()